rm(list = ls())

# This creates the first part of figures and tables in Appendix

# Load Packages
library(tidyverse) # 1.3.1
library(stargazer) # 5.2.2
library(scales) # 1.1.1
library(modelsummary) # 0.7.0
library(vtable) # 1.3.1
library(gtsummary) # 1.4.1
library(gridExtra) # 2.3
library(jtools) # 2.1.3
library(sandwich) # 3.0-1


SEs <- function(data){
  empty_data <- c()
  data1 <- vcovHC(data, type = "HC1")
  empty_data <- sqrt(diag(data1))
  return(empty_data)
}

### Load data
data <- read.csv("perspective_taking_merged.csv")

## Filter out those who refused informed consent 
# + failed attention checks
data %>% 
  filter(informed_consent == 1 & attention_check1 == 3) %>%
  mutate_all(na_if,"") %>%
  filter(age >= 18) %>%
  filter(!is.na(attention_check1)) -> data


######################################   
########## Figure A4 #################
######################################


###### Plot Including Low Income Group
data %>%
  mutate(order = recode(treatment, ## Provide an ordering label
                        "Control Undoc" = 1,
                        "Treat Undoc" = 2,
                        "Cancer" = 3,
                        "Immigrant" = 4,
                        "Control Low Inc" = 5,
                        "Treat Low Inc" = 6)) %>%
  group_by(pid_3level, treatment, treated, order) %>%
  mutate(treatment = factor(treatment, # Arrange in the right order
                            levels = c("Control Undoc",
                                       "Treat Undoc",
                                       "Cancer",
                                       "Immigrant",
                                       "Control Low Inc",
                                       "Treat Low Inc"))) %>%
  summarise(n = n(),
            mean = mean(response_deserving, na.rm = T),
            se = sd(response_deserving, na.rm =T)/sqrt(n())) %>%
  filter(!is.na(pid_3level),
         !is.na(treatment)) %>%
  ggplot(aes(x = reorder(treatment, -order), y = mean,
             ymax = (mean + 2*se), 
             ymin = (mean - 2*se),
             color = as.factor(pid_3level), 
             shape = as.factor(pid_3level)) ) +
  geom_pointrange(position = position_dodge2(w = 0.5)) +
  coord_flip() +
  theme_bw() +
  scale_y_continuous(labels = scales::percent) +
  scale_x_discrete(labels = c("Control Undoc" = "Control\nUndocumented",
                              "Treat Undoc" = "Perspective Taking\nUndocumented\n+COVID-19",
                              "Cancer" = "Perspective Taking\nUndocumented\n+Cancer",
                              "Immigrant" = "Perspective Taking\nImmigrant\n+COVID-19",
                              "Control Low Inc" = "Control\nLow Income",
                              "Treat Low Inc" = "Perspective Taking\nLow Income\n+COVID-19") ) +
  scale_shape_discrete(name = "",
                       labels = c("Democrats",
                                  "Independents",
                                  "Republicans")) +
  scale_color_manual(name = "",
                     labels = c("Democrats",
                                "Independents",
                                "Republicans"),
                     values = c("skyblue3",
                                "purple3",
                                "red3")) +
  theme(legend.position = "top",
        axis.text.x = element_text(size = rel(1.45)),
        axis.text.y = element_text(size = rel(1.45))) +
  labs(x = "", y = "", title = "")
ggsave("appendix_figure_a4_final.jpeg")



######################################
######### Table A6 and A7 ############
######################################

# Study 1 - Immigrant groups
study1_immig <- data %>%
  filter(study == 1,
         !treatment %in% c("Control Low Inc",
                           "Treat Low Inc") ) %>%
  droplevels()

# Study 1 - Low Income Groups
study1_lowinc <- data %>%
  filter(study == 1,
         treatment %in% c("Control Low Inc",
                           "Treat Low Inc") ) %>%
  droplevels()

# Study 2
study2 <- data %>%
  filter(study == 2) %>%
  mutate(treatment = factor(treatment,
                            levels = c("Control Undoc",
                                       "Treat Undoc",
                                       "Cancer",
                                       "Immigrant")))


### Create Table A6
# Get regressions + SEs
deserving_study1_immig <- lm(response_deserving ~ treatment, study1_immig)
deserving_study1_immig_se <- as.data.frame(SEs(deserving_study1_immig))

deserving_study1_lowinc <- lm(response_deserving ~ treatment, study1_lowinc)
deserving_study1_lowinc_se <- as.data.frame(SEs(deserving_study1_lowinc))

deserving_study2 <- lm(response_deserving ~ treatment, study2)
deserving_study2_se <- as.data.frame(SEs(deserving_study2))

# Make Table
stargazer(deserving_study1_immig,
          deserving_study2,
          deserving_study1_lowinc,  
          type = "text", align = T,
          font.size = "footnotesize",
          digits = 3, 
          no.space = T,
          initial.zero = F,
          keep = c("Constant", "treatmentTreat Undoc",
                   "treatmentCancer", "treatmentImmigrant",
                   "treatmentTreat Low Inc"),
          se = list(deserving_study1_immig_se[,1], 
                    deserving_study2_se[,1],
                    deserving_study1_lowinc_se[,1]),
          column.labels = c("Study 1", "Study 2", "Study 1"),
          model.numbers = F,
          covariate.labels = c("Undocumented + COVID-19",
                               "Undocumented + Cancer",
                               "Immigrant + COVID-19",
                               "Low Income + COVID-19",
                               "Constant"),
          dep.var.labels = "Deserving of Gov Healthcare")


### Create Table A7

# Study 1 - Undoc Immigrants as targed group
deserving_study1_immig_adj <- lm(response_deserving ~ treatment
                                 + age + gender + income + educ + ethnicity, study1_immig)
deserving_study1_immig_adj_se <- as.data.frame(SEs(deserving_study1_immig_adj))

# Study 1 - Low Income individuals as targeted group
deserving_study1_lowinc_adj <- lm(response_deserving ~ treatment
                              + age + gender + income + educ + ethnicity, study1_lowinc)
deserving_study1_lowinc_adj_se <- as.data.frame(SEs(deserving_study1_lowinc_adj))

# Study 2 - Immigrants + Undocumented immigrants as targeted group
deserving_study2_adj <- lm(response_deserving ~ treatment
                           + age + gender + income + educ + ethnicity, study2)
deserving_study2_adj_se <- as.data.frame(SEs(deserving_study2_adj))

### Create Table
stargazer(deserving_study1_immig_adj,
          deserving_study2_adj,
          deserving_study1_lowinc_adj,  
          type = "text", align = T,
          font.size = "footnotesize",
          digits = 3, 
          no.space = T,
          initial.zero = F,
          keep = c("Constant", "treatmentTreat Undoc",
                   "treatmentCancer", "treatmentImmigrant",
                   "treatmentTreat Low Inc"),
          se = list(deserving_study1_immig_adj_se[,1], 
                    deserving_study2_adj_se[,1],
                    deserving_study1_lowinc_adj_se[,1]),
          column.labels = c("Study 1", "Study 2", "Study 1"),
          model.numbers = F,
          covariate.labels = c("Undocumented + COVID-19",
                               "Undocumented + Cancer",
                               "Immigrant + COVID-19",
                               "Low Income + COVID-19",
                               "Constant"),
          dep.var.labels = "Deserving of Gov Healthcare")


######################################
######### Table A8 and A9 ############
######################################


## Make Plots Across Party ID 
data %>%
  filter(treatment %in% c("Control Undoc", 
                          "Treat Undoc", # Focus on groups with undocumented/immigrants
                          "Cancer", 
                          "Immigrant")) %>%
  mutate(treatment = factor(treatment, # Arrange in a better order
                             levels = c("Control Undoc",
                                        "Treat Undoc",
                                        "Cancer",
                                        "Immigrant"))) -> study_merge_table_pid
## Subset by Party ID
# Republicans
reps <- study_merge_table_pid %>%
  filter(pid_3level == 3,
         !treatment %in% c("Control Low Inc",
                            "Treat Low Inc") )
# Democrats
dems <- study_merge_table_pid %>%
  filter(pid_3level == 1,
         !treatment %in% c("Control Low Inc",
                            "Treat Low Inc"))
# Independents
indep <- study_merge_table_pid %>%
  filter(pid_3level == 2,
         !treatment %in% c("Control Low Inc",
                            "Treat Low Inc"))


######################################
############ TABLE A8 ################
########### Unadjusted ###############
######################################

# Produce table with unadjusted results
# viewer
modelsummary(list("Pooled\nUnadjusted" = lm(response_deserving ~ treatment, study_merge_table_pid),
                  "Democrats\nUnadjusted" = lm(response_deserving ~ treatment, dems),
                  "Independents\nUnadjusted" = lm(response_deserving ~ treatment, indep),
                  "Republicans\nUnadjusted" = lm(response_deserving ~ treatment, reps)),
             coef_rename = c("treatmentTreat Undoc"= "Perspective Taking\nCOVID-19 + Undoc", ### Pooled Study 1 and Study 2
                             "treatmentCancer" = "Perspective Taking\nCancer + Undocumented", ## Study 2
                             "treatmentImmigrant" = "Perspective Taking\nCOVID-19 + Immigrant"), ## Study 2
             vcov = "HC1", stars = T,
             leading_zero = F)

# using stargazer
# Pooled unadjusted
deserving_pooled_unadj <- lm(response_deserving ~ treatment, study_merge_table_pid)
deserving_pooled_unadj_se <- as.data.frame(SEs(deserving_pooled_unadj))

# Republican unadjusted
deserving_rep_unadj <- lm(response_deserving ~ treatment, reps)
deserving_rep_unadj_se <- as.data.frame(SEs(deserving_rep_unadj))

# Independent unadjusted
deserving_indep_unadj <- lm(response_deserving ~ treatment, indep)
deserving_indep_unadj_se <- as.data.frame(SEs(deserving_indep_unadj))

# Democrats unadjusted
deserving_dem_unadj <- lm(response_deserving ~ treatment, dems)
deserving_dem_unadj_se <- as.data.frame(SEs(deserving_dem_unadj))

# Make Table
stargazer(deserving_pooled_unadj,
          deserving_dem_unadj,
          deserving_indep_unadj, 
          deserving_rep_unadj,
          type = "text", align = T,
          font.size = "footnotesize",
          digits = 3, 
          no.space = T,
          initial.zero = F,
          keep = c("Constant", "treatmentTreat Undoc",
                   "treatmentCancer", "treatmentImmigrant",
                   "treatmentTreat Low Inc"),
          se = list(deserving_pooled_unadj_se[,1],
                    deserving_dem_unadj_se[,1], 
                    deserving_indep_unadj_se[,1],
                    deserving_rep_unadj_se[,1]),
          column.labels = c("Pooled", "Democrats", 
                            "Independents", "Republicans"),
          model.numbers = F,
          covariate.labels = c("Undocumented + COVID-19",
                               "Undocumented + Cancer",
                               "Immigrant + COVID-19",
                               "Constant"),
          dep.var.labels = "Undocumented Immigrants as Deserving of Gov Healthcare")

######################################
############ TABLE A9 ################
############ Adjusted ################
######################################

# viewer
modelsummary(list("Pooled\nAdjusted" = lm(response_deserving ~ treatment
                                            + age + gender + income + educ + ethnicity, 
                                            study_merge_table_pid),
                  "Democrats\nAdjusted" = lm(response_deserving ~ treatment
                                               + age + gender + income + educ + ethnicity,
                                               dems),
                  "Independents\nAdjusted" = lm(response_deserving ~ treatment
                                                  + age + gender + income + educ + ethnicity,
                                                  indep),
                  "Republicans\nAdjusted" = lm(response_deserving ~ treatment
                                                 + age + gender + income + educ + ethnicity, 
                                                 reps)),
             coef_rename = c("treatmentTreat Undoc"= "Perspective Taking\nCOVID-19\n+ Undocumented", ### Pooled Study 1 and Study 2
                             "treatmentCancer" = "Perspective Taking\nCancer\n+ Undocumented", ## Study 2
                             "treatmentImmigrant" = "Perspective Taking\nCOVID-19\n+ Immigrant"), ## Study 2
             vcov = "HC1", stars = T,
             leading_zero = F)


### Produce Table adjusted
## Pooled
deserving_pooled_adj <- lm(response_deserving ~ treatment 
                       + age + educ + income + gender + ethnicity, study_merge_table_pid)
deserving_pooled_adj_se <- as.data.frame(SEs(deserving_pooled_adj))

# Republicans
deserving_rep_adj <- lm(response_deserving ~ treatment 
                    + age + educ + income + gender + ethnicity, reps)
deserving_rep_adj_se <- as.data.frame(SEs(deserving_rep_adj))

# Independents
deserving_indep_adj <- lm(response_deserving ~ treatment 
                      + age + educ + income + gender + ethnicity, indep)
deserving_indep_adj_se <- as.data.frame(SEs(deserving_indep_adj))

# Democrats
deserving_dem_adj <- lm(response_deserving ~ treatment 
                    + age + educ + income + gender + ethnicity, dems)
deserving_dem_adj_se <- as.data.frame(SEs(deserving_dem_adj))

# Export
stargazer(deserving_pooled_adj,
          deserving_dem_adj, 
          deserving_indep_adj, 
          deserving_rep_adj,
          type = "text", align = T,
          font.size = "footnotesize",
          digits = 3, 
          no.space = T,
          initial.zero = F,
          keep = c("Constant", 
                   "treatmentTreat Undoc",
                   "treatmentCancer", 
                   "treatmentImmigrant"),
          se = list(deserving_pooled_adj_se[,1], 
                    deserving_dem_adj_se[,1], 
                    deserving_indep_adj_se[,1], 
                    deserving_rep_adj_se[,1]),
          column.labels = c("Pooled", "Democrats", 
                            "Independents", "Republicans"),
          model.numbers = F,
          covariate.labels = c("Undocumented + COVID-19",
                               "Undocumented + Cancer",
                               "Immigrant + COVID-19",
                               "Constant"),
          dep.var.labels = "Undocumented Immigrants as Deserving of Gov Healthcare")

######################################
########## Figure A5 #################
######################################

## Create figure across favorability of undocumented immigrants
data %>%
  filter(treatment %in% c("Control Undoc", ## Focus on the undocumented groups
                          "Treat Undoc", 
                          "Cancer")) %>%
  mutate(treatment = factor(treatment, # Arrange in a better order
                              levels = c("Control Undoc",
                                         "Treat Undoc",
                                         "Cancer"))) %>%
  group_by(undoc_immigrants, treatment) %>%
  summarise(sample_count = n(),
            mean = mean(response_deserving, na.rm = T),
            se = sd(response_deserving, na.rm = T) / sqrt(n())) %>%
  filter(!is.na(undoc_immigrants),
         !is.na(treatment)) %>%
  ggplot(aes(x = treatment, y = mean,
             ymin = mean - (2*se), 
             ymax = mean + (2*se)),
         shape = as.factor(treatment), 
         color = as.factor(treatment)) +
  geom_pointrange(position = position_dodge2(w = 0.5)) +
  facet_wrap(undoc_immigrants ~.,
             nrow = 5,
             scales = "free",
             labeller = as_labeller(c(`0` = "0 - Very unfavorable",
                                      `1` = "1",
                                      `2` = "2",
                                      `3` = "3",
                                      `4` = "4 - Very favorable")) ) +
  labs(x = "", 
       y = "Mean Deservingness Support") +
  scale_y_continuous(labels = scales::percent) +
  scale_x_discrete(labels = c("Control", "Perspective Taking\n+COVID-19", "Perspective Taking\n+Cancer")) +
  theme_bw() + 
  theme(legend.position = "top") 
ggsave("appendix_figure_a5_final.jpeg")

#########################################
############ Policy #####################
#########################################

######################################
########## Figure A6 #################
######################################

# Create figure with policy results in study 1 
data %>%
  filter(study == 1) %>% # Filter for study 1
  mutate(treated = case_when(treatment == "Control Undoc" |
                             treatment == "Control Low Inc" ~ "Control", # Label control and treatment groups
                             TRUE ~ "Perspective Taking"),
         group = case_when(treatment == "Control Low Inc" |
                             treatment == "Treat Low Inc" ~ "Low Inc", # Label low income/undocumented groups
                           TRUE ~ "Undocumented")) %>%
  group_by(treatment, treated, group) %>%
  summarise(n = n(),
            mean = mean(response_policy, na.rm = T),
            se = sd(response_policy, na.rm = T) / sqrt(n())) %>%
  filter(!is.na(treatment)) %>%
  ggplot(aes(x = treatment, y = mean,
             ymax = mean + (2*se), 
             ymin = mean - (2*se),
             shape = treated, color = treated)) +
  geom_pointrange(position = position_dodge(w = 0.75)) +
  labs(x = "", y = "Mean Policy Support") +
  scale_y_continuous(labels = scales::percent) +
  scale_x_discrete(labels=c('Control', 'Perspective-Taking')) +
  theme_bw() +
  theme(legend.position = "top",
        axis.text.x=element_text(size=rel(1.5)),
        axis.text.y = element_text(size = rel(1.5))) +
  coord_flip() +
  scale_color_manual(name = "",
                     labels = c("Control", "Perspective-Taking"),
                     values = c("dark blue", "red3")) +
  scale_shape_discrete(name = "",
                       labels = c("Control", "Perspective-Taking")) +
  facet_wrap(. ~ group, 
             strip.position = "top",
             dir = "v",
             scales = "free_y",
             labeller = as_labeller(c(`Low Inc` = "Low Income",
                                      `Undocumented` = "Undocumented Immigrant")) ) 
ggsave("appendix_figure_a6_final.jpeg")


######################################
####### Policy Results ##################
######################################

# Select only the undoc groups in study 1
policy_table_undoc <- data %>%
  filter(study == 1 & ## Select only Study 1
           treatment %in% c("Control Undoc", ### Select only groups focused on Undoc
                            "Treat Undoc")) %>%
  mutate(treatment = factor(treatment, ## Set Control as reference
                            levels = c("Control Undoc",
                                       "Treat Undoc")),
         treated = ifelse(treatment == "Control Undoc",
                          "Control", "Perspective Taking")) 

# Select only the low inc groups in study 1
policy_table_lowinc <- data %>%
  filter(study == 1 & ## Select only Study 1
           !treatment %in% c("Control Undoc", ### Select only groups focused on Low Inc
                            "Treat Undoc")) %>%
  mutate(treatment = factor(treatment, ## Set Control as reference
                            levels = c("Control Low Inc",
                                       "Treat Low Inc")),
         treated = ifelse(treatment == "Control Low Inc",
                          "Control", "Perspective Taking")) 

######################################
####### Table A10 ##################
######################################

# Policy Undoc Unadjusted
policy_undoc_unadj <- lm(response_policy ~ treatment, policy_table_undoc)
policy_undoc_unadj_se <- as.data.frame(SEs(policy_undoc_unadj))

# Policy Undoc Adjusted
policy_undoc_adj <- lm(response_policy ~ treatment
                       + age + educ + income + gender + ethnicity, policy_table_undoc)
policy_undoc_adj_se <- as.data.frame(SEs(policy_undoc_adj))

# Policy Low Inc Unadjusted
policy_lowinc_unadj <- lm(response_policy ~ treatment, policy_table_lowinc)
policy_lowinc_unadj_se <- as.data.frame(SEs(policy_lowinc_unadj))


# Policy Low Inc Adjusted
policy_lowinc_adj <- lm(response_policy ~ treatment
                       + age + educ + income + gender + ethnicity, policy_table_lowinc)
policy_lowinc_adj_se <- as.data.frame(SEs(policy_lowinc_adj))

# Create Table
modelsummary(list(
  "Undoc\nUnadjusted" = policy_undoc_unadj,
  "Undoc\nAdjusted" = policy_undoc_adj,
  "Low Income\nUnadjusted" = policy_lowinc_unadj,
  "Low Income\nAdjusted" = policy_lowinc_adj),
  stars = T, 
  vcov = "HC1",
  coef_rename = c("treatmentTreat Undoc" = "Perspective Taking\nUndocumented",
                  "treatmentTreat Low Inc" = "Perspective Taking\nLow Income"),
  out = "default")


# Create Table
stargazer(policy_undoc_unadj,
          policy_undoc_adj,
          policy_lowinc_unadj,
          policy_lowinc_adj,
          type = "text", align = T,
          font.size = "footnotesize",
          digits = 3, 
          no.space = T,
          initial.zero = F,
          keep = c("Constant",
                   "treatmentTreat Undoc",
                   "treatmentTreat Low Inc"),
          se = list(policy_undoc_unadj_se[,1], 
                    policy_undoc_adj_se[,1], 
                    policy_lowinc_unadj_se[,1], 
                    policy_lowinc_adj_se[,1]),
          column.labels = c("Undoc Unadjusted", "Undoc Adjusted", 
                            "Low Income Unadjusted", "Low Income Adjusted"),
          model.numbers = F,
          covariate.labels = c("Undocumented + COVID-19",
                               "Low Income + COVID-19",
                               "Constant"),
          dep.var.labels = c("Policy Preferences"))




######################################
########## Tables A11 and A12 ########
######################################

# Prepare data for table
policy_table <- data %>%
  filter(study == 1 & ## Select only Study 1
           treatment %in% c("Control Undoc", ### Select only Treatment focused on Undoc
                            "Treat Undoc")) %>%
  mutate(treatment = factor(treatment, ## Set Control as reference
                            levels = c("Control Undoc",
                                       "Treat Undoc")),
         treated = ifelse(treatment == "Control Undoc",
                          "Control", "Perspective Taking")) 

# Democrats
dems <- policy_table %>%
  filter(pid_3level == 1)

# Independents
indep <- policy_table %>%
  filter(pid_3level == 2)

# Republicans
reps <- policy_table %>%
  filter(pid_3level == 3)


######################################
########## Table A11 ##################
######################################


# Create Table with unadjusted results
# Pooled
policy_pooled_unadj <- lm(response_policy ~ treatment, policy_table)
policy_pooled_unadj_se <- as.data.frame(SEs(policy_pooled_unadj))

# Democrats
policy_dems_unadj <- lm(response_policy ~ treatment, dems)
policy_dems_unadj_se <- as.data.frame(SEs(policy_dems_unadj))

# Independents
policy_indep_unadj <- lm(response_policy ~ treatment, indep)
policy_indep_unadj_se <- as.data.frame(SEs(policy_indep_unadj))

# Republicans
policy_reps_unadj <- lm(response_policy ~ treatment, reps)
policy_reps_unadj_se <- as.data.frame(SEs(policy_reps_unadj))

# Create Table
modelsummary(list(
  "Pooled\nUnadjusted" = policy_pooled_unadj,
  "Democrats\nUnadjusted" = policy_dems_unadj,
  "Independents\nUnadjusted" = policy_indep_unadj,
  "Republicans\nUnadjusted" = policy_reps_unadj),
  stars = T, 
  vcov = "HC1",
  coef_rename = c("treatmentTreat Undoc" = "Perspective Taking\nUndocumented"),
  out = "default")


# Create Table
stargazer(policy_pooled_unadj,
          policy_dems_unadj,
          policy_indep_unadj,
          policy_reps_unadj,
          type = "text", align = T,
          font.size = "footnotesize",
          digits = 3, 
          no.space = T,
          initial.zero = F,
          keep = c("Constant",
                   "treatmentTreat Undoc"),
          se = list(policy_pooled_unadj_se[,1], 
                    policy_dems_unadj_se[,1], 
                    policy_indep_unadj_se[,1], 
                    policy_reps_unadj_se[,1]),
          column.labels = c("Pooled", "Democrats", 
                            "Independents", "Republicans"),
          model.numbers = F,
          covariate.labels = c("Undocumented + COVID-19",
                               "Constant"),
          dep.var.labels = c("Policy Preferences"))




######################################
########## Table A12 ##################
######################################

# Create Table with adjusted results
# Pooled
policy_pooled_adj <- lm(response_policy ~ treatment 
                          + age + educ + income + gender + ethnicity, policy_table)
policy_pooled_adj_se <- as.data.frame(SEs(policy_pooled_adj))

# Democrats
policy_dems_adj <- lm(response_policy ~ treatment 
                        + age + educ + income + gender + ethnicity, dems)
policy_dems_adj_se <- as.data.frame(SEs(policy_dems_adj))

# Independents
policy_indep_adj <- lm(response_policy ~ treatment 
                         + age + educ + income + gender + ethnicity, indep)
policy_indep_adj_se <- as.data.frame(SEs(policy_indep_adj))

# Republicans
policy_reps_adj <- lm(response_policy ~ treatment 
                        + age + educ + income + gender + ethnicity, reps)
policy_reps_adj_se <- as.data.frame(SEs(policy_reps_adj))

# Create Table
modelsummary(list(
  "Pooled\nAdjusted" = policy_pooled_adj,
  "Democrats\nAdjusted" = policy_dems_adj,
  "Independents\nAdjusted" = policy_indep_adj,
  "Republicans\nAdjusted" = policy_reps_adj),
  stars = T, 
  vcov = "HC1",
  coef_rename = c("treatmentTreat Undoc" = "Perspective Taking\nUndocumented"),
  out = "default")


# Create Table
stargazer(policy_pooled_adj,
          policy_dems_adj,
          policy_indep_adj,
          policy_reps_adj,
          type = "text", align = T,
          font.size = "footnotesize",
          digits = 3, 
          no.space = T,
          initial.zero = F,
          keep = c("Constant",
                   "treatmentTreat Undoc"),
          se = list(policy_pooled_adj_se[,1], 
                    policy_dems_adj_se[,1], 
                    policy_indep_adj_se[,1], 
                    policy_reps_adj_se[,1]),
          column.labels = c("Pooled", 
                            "Democrats", 
                            "Independents", 
                            "Republicans"),
          model.numbers = F,
          covariate.labels = c("Undocumented + COVID-19",
                               "Constant"),
          dep.var.labels = c("Policy Preferences"))




######################################
########### Study 2 Policy ###########
######################################

policy_table_exp2 <- data %>%
  filter(study == 2) %>% ### Select only Study 2
  mutate(treatment_s2 = factor(treatment_s2,  ## Select only the Second Policy Experiment 
                               levels = c("Health",
                                          "Health & Welfare",
                                          "Welfare")),
         treated = factor(case_when(treatment == "Control Undoc" ~ "Control", ## Label the Control 
                             treatment == "Treat Undoc" ~ "Treatment Undoc", ## Label Treat + Undoc + COVID-19
                             treatment == "Cancer" ~ "Cancer", ## Label Treat + Undoc + Cancer
                             treatment == "Immigrant" ~ "Immigrant"),
                          levels = c(levels = c("Control",
                               "Treatment Undoc",
                               "Cancer",
                               "Immigrant"))) )


######################################
####### Tables A13 and A14 ###########
######################################


# Health 
health <- policy_table_exp2 %>%
  filter(treatment_s2 == "Health")

# Welfare
welfare <- policy_table_exp2 %>%
  filter(treatment_s2 == "Welfare")

# Health and Welfare
health_welfare <- policy_table_exp2 %>%
  filter(treatment_s2 == "Health & Welfare")

######################################
########## Table A13 ##################
######################################


# Health unadjusted
policy_table_health_unadj <- lm(response_policy ~ treated, health) 
policy_table_health_unadj_se <- as.data.frame(SEs(policy_table_health_unadj))

# Welfare unadjusted
policy_table_welfare_unadj <- lm(response_policy ~ treated, welfare)
policy_table_welfare_unadj_se <- as.data.frame(SEs(policy_table_welfare_unadj))

# Health and Welfare unadjusted
policy_table_health_welfare_unadj <- lm(response_policy ~ treated, health_welfare) 
policy_table_health_welfare_unadj_se <- as.data.frame(SEs(policy_table_health_welfare_unadj))

## Create Table Unadjusted
modelsummary(list(
  "Health\nUnadjusted" = policy_table_health_unadj,
  "Health & Welfare\nUnadjusted" = policy_table_health_welfare_unadj,
  "Welfare\nUnadjusted" = policy_table_welfare_unadj),
  stars = T, 
  vcov = "HC1",
  coef_rename = c("treatedTreatment Undoc"= "Perspective Taking\n+COVID-19 + Undoc", ### Pooled Study 1 and Study 2
                  "treatedCancer" = "Perspective Taking\n+Cancer + Undocumented", ## Study 2
                  "treatedImmigrant" = "Perspective Taking\n+COVID-19 + Immigrant"),
  out = "default")

## Create Table Adjusted
stargazer(policy_table_health_unadj,
          policy_table_health_welfare_unadj,
          policy_table_welfare_unadj,
          type = "text", align = T,
          font.size = "footnotesize",
          digits = 3, 
          no.space = T,
          initial.zero = F,
          keep = c("Constant",
                   "treatedTreatment Undoc",
                   "treatedCancer",
                   "treatedImmigrant"),
          se = list(policy_table_health_unadj_se[,1], 
                    policy_table_health_welfare_unadj_se[,1], 
                    policy_table_welfare_unadj_se[,1]),
          column.labels = c("Health", 
                            "Health + Welfare", 
                            "Welfare"),
          model.numbers = F,
          covariate.labels = c("Undocumented + COVID-19",
                               "Undocumented + Cancer",
                               "Immigrant + COVID-19",
                               "Constant"),
          dep.var.labels = c("Policy Preferences"))



######################################
########## Table A14 ##################
######################################


# Health adjusted
policy_table_health_adj <- lm(response_policy ~ treated
                               + age + educ + income + gender + ethnicity, health) 
policy_table_health_adj_se <- as.data.frame(SEs(policy_table_health_adj))

# Welfare adjusted
policy_table_welfare_adj <- lm(response_policy ~ treated
                                + age + educ + income + gender + ethnicity, welfare)
policy_table_welfare_adj_se <- as.data.frame(SEs(policy_table_welfare_adj))

# Health and Welfare adjusted
policy_table_health_welfare_adj <- lm(response_policy ~ treated
                                       + age + educ + income + gender + ethnicity, health_welfare) 
policy_table_health_welfare_adj_se <- as.data.frame(SEs(policy_table_health_welfare_adj))


modelsummary(list(
  "Health\nAdjusted" = policy_table_health_adj,
  "Health & Welfare\nAdjusted" = policy_table_health_welfare_adj,
  "Welfare\nAdjusted" = policy_table_welfare_adj),
  stars = T, 
  vcov = "HC1",
  coef_rename = c("treatedTreatment Undoc"= "Perspective Taking\n+COVID-19 + Undoc", ### Pooled Study 1 and Study 2
                  "treatedCancer" = "Perspective Taking\n+Cancer + Undocumented", ## Study 2
                  "treatedImmigrant" = "Perspective Taking\n+COVID-19 + Immigrant"),
  out = "default")


## Create Table Adjusted
stargazer(policy_table_health_adj,
          policy_table_health_welfare_adj,
          policy_table_welfare_adj,
          type = "text", align = T,
          font.size = "footnotesize",
          digits = 3, 
          no.space = T,
          initial.zero = F,
          keep = c("Constant",
                   "treatedTreatment Undoc",
                   "treatedCancer",
                   "treatedImmigrant"),
          se = list(policy_table_health_adj_se[,1], 
                    policy_table_health_welfare_adj_se[,1], 
                    policy_table_welfare_adj_se[,1]),
          column.labels = c("Health", 
                            "Health + Welfare", 
                            "Welfare"),
          model.numbers = F,
          covariate.labels = c("Undocumented + COVID-19",
                               "Undocumented + Cancer",
                               "Immigrant + COVID-19",
                               "Constant"),
          dep.var.labels = c("Policy Preferences"))


#####################################   
######## Figure A7 ################
#####################################   

data %>%
  filter(study == 2) %>% # Select Study 2
  mutate(treatment = factor(treatment, # Reorder treatment in better order
                            levels = c("Control Undoc",
                                       "Treat Undoc",
                                       "Cancer",
                                       "Immigrant")),
         treatment_s2 = factor(treatment_s2, # Reorder second survey experiment
                               levels = c("Welfare",
                                          "Health",
                                          "Health & Welfare")),
         treated = ifelse(treatment == "Control Undoc", 
                          "Control", "Perspective Taking")) %>%
  group_by(treatment_s2, treated, treatment) %>%
  summarise(n = n(),
            mean = mean(response_policy, na.rm = T),
            se = sd(response_policy, na.rm = T) / sqrt(n())) %>%
    filter(!is.na(treated),
           !is.na(treatment_s2)) %>%
  ggplot(aes(x = treatment_s2, 
             y = mean, 
             ymax = mean + (2*se),
             ymin = mean - (2*se),
             shape = as.factor(treated),
             color = as.factor(treated))) +
  geom_pointrange(position = position_dodge2(w = .5)) +
  theme_bw() +
  facet_wrap(.~ treatment,
             scales = "free_x",
             labeller = as_labeller(c(`Control Undoc` = "Control",
                                      `Treat Undoc` = "Perspective-Taking \nCOVID-19",
                                      `Cancer` = "Perspective-Taking \nCancer",
                                      `Immigrant` = "Perspective-Taking \nImmigrant") ) ) +
  scale_y_continuous(labels = scales::percent) +
  scale_color_manual(name = "",
                     labels = c("Control",
                                "Perspective-Taking"),
                     values = c("dark blue",
                                "red3")) +
  scale_shape_discrete(name = "",
                       labels = c("Control",
                                  "Perspective-Taking")) +
  theme(legend.position = "top",
        axis.text.x = element_text(size = rel(1.25)),
        axis.text.y = element_text(size = rel(1.25))) +
  scale_x_discrete(labels = c("Welfare" = "Welfare",
                              "Health" = "Health",
                              "Health & Welfare" = "Health and Welfare") ) +
  labs(x = "", y = "Mean Policy Support", title = "")
ggsave("appendix_figure_a7_final.jpeg")


#####################################   
######## Figure A8 ################
#####################################

# Label COVID-19 Lived Experience
data %>%
  mutate(covid_ever = ifelse(covid_contracted == 1 |
                               covid_self == 1, "Yes", "No")) -> data

# Create Figure of COVID-19 Lived Experience across treatment groups
data %>% 
  mutate(treatment = factor(treatment, # Reorder treatment variable
                            levels = c("Control Undoc",
                                       "Treat Undoc",
                                       "Cancer",
                                       "Immigrant",
                                       "Control Low Inc",
                                       "Treat Low Inc")),
         order = recode(treatment, # Create ordering indicator
                        "Control Undoc" = 1,
                        "Treat Undoc" = 2,
                        "Cancer" = 3,
                        "Immigrant" = 4,
                        "Control Low Inc" = 5,
                        "Treat Low Inc" = 6)) %>% 
  group_by(treatment, covid_ever, treated, order) %>%
  summarise(n = n(),
            mean = mean(response_deserving, na.rm = T),
            se = sd(response_deserving, na.rm = T) / sqrt(n())) %>%
  filter(!is.na(treatment),
         !is.na(covid_ever)) %>%
  ggplot(aes(x = reorder(treatment, -order), y = mean,
             ymax = (mean + 2*se), 
             ymin = (mean - 2*se),
             shape = as.factor(covid_ever), 
             color = as.factor(covid_ever)) ) +
  geom_pointrange(position = position_dodge2(w = 0.5)) +
  coord_flip() +
  theme_bw() +
  scale_y_continuous(labels = scales::percent) +
  scale_shape_discrete(name = "COVID-19 \nPersonal Experience",
                       labels = c("No",
                                  "Yes")) +
  scale_color_manual(name = "COVID-19 \nPersonal Experience",
                       labels = c("No",
                                  "Yes"),
                       values = c("dark blue",
                                  "red3")) +
  scale_x_discrete(labels = c("Control Undoc" = "Control\n+Undocumented",
                              "Treat Undoc" = "Perspective Taking\n+Undocumented\n+COVID-19",
                              "Cancer" = "Perspective Taking\n+Undocumented\n+Cancer",
                              "Immigrant" = "Perspective Taking\n+Immigrant\n+COVID-19",
                              "Control Low Inc" = "Control\n+Low Income",
                              "Treat Low Inc" = "Perspective Taking\n+Low Income\n+COVID-19") ) + 
  theme(legend.position = "top",
        axis.text.x = element_text(size = rel(1.35)),
        axis.text.y = element_text(size = rel(1.15))) +
  labs(x = "", y = "Mean Deservingness Support", title = "")
ggsave("appendix_figure_a8_final.jpeg")


#####################################   
######## Figure A9 ################
#####################################

# Create figure of COVID-19 Lived Experience across party-ID
data %>% 
  filter(treatment %in% c("Control Undoc", # subset among groups with undocumented
                          "Treat Undoc")) %>%
  mutate(treated = factor(treated,  # Create indicator for treated/control
                          levels = c("Treated",
                                     "Control"))) %>%
  group_by(covid_ever, treated, pid_3level) %>%
  summarise(n = n(),
            mean = mean(response_deserving, na.rm = T),
            se = sd(response_deserving, na.rm = T) / sqrt(n())) %>%
  filter(!is.na(covid_ever),
         !is.na(pid_3level)) %>%
  ggplot(aes(x = treated, y = mean,
             ymax = (mean + 2*se), 
             ymin = (mean - 2*se),
             shape = as.factor(covid_ever), 
             color = as.factor(covid_ever)) ) +
  geom_pointrange(position = position_dodge2(w = 0.5)) +
  theme_bw() +
  scale_y_continuous(labels = scales::percent) +
  facet_wrap(.~ pid_3level,
             dir = "v",
             labeller = as_labeller(c(`1` = "Democrats",
                                      `2` = "Independents",
                                      `3` = "Republicans"))) +
  scale_x_discrete(labels = c("Control" = "Control",
                              "Treated" = "Perspective Taking") ) +
  coord_flip() +
  theme_bw() +
  scale_color_manual(name = "COVID-19\nLived Experience",
                     labels = c("No",
                                "Yes"),
                     values = c("dark blue",
                                "red3")) +
  scale_shape_discrete(name = "COVID-19\nLived Experience", 
                       labels = c("No",
                                  "Yes")) +
  theme(legend.position = "top") +
  labs(x = "", y = "Mean Deservingness Support", title = "")
ggsave("appendix_figure_a9_final.jpeg")



#####################################   
######## Figure A10 ################
#####################################

# Create Figure for COVID-19 personal + family Lived Experience
data %>% 
  mutate(covid_ever_fam = case_when(
    covid_family == 1 | covid_ever == 1 ~ "Yes", TRUE ~ "No")) %>% # create indicator
  filter(treatment %in% c("Control Undoc", # subset among groups with undocumented
                          "Treat Undoc")) %>%
  mutate(treated = factor(treated, # Create indicator for treated/control
                          levels = c("Treated",
                                     "Control"))) %>%
  group_by(covid_ever_fam, treated) %>%
  summarise(n = n(),
            mean = mean(response_deserving, na.rm = T),
            se = sd(response_deserving, na.rm = T) / sqrt(n())) %>%
  filter(!is.na(covid_ever_fam)) %>%
  ggplot(aes(x = treated, y = mean,
             ymax = (mean + 2*se), 
             ymin = (mean - 2*se),
             shape = as.factor(covid_ever_fam), 
             color = as.factor(covid_ever_fam)) ) +
  geom_pointrange(position = position_dodge2(w = 0.5)) +
  coord_flip() +
  theme_bw() +
  scale_y_continuous(labels = scales::percent) +
  scale_x_discrete(labels = c("Control" = "Control",
                              "Treated" = "Perspective Taking") ) +
  scale_color_manual(name = "COVID-19 \nPersonal + Family \nExperience",
                     labels = c("No",
                                "Yes"),
                     values = c("dark blue",
                                "red3")) +
  scale_shape_discrete(name = "COVID-19 \nPersonal + Family \nExperience",
                       labels = c("No",
                                  "Yes")) +
  theme(legend.position = "top") +
  labs(x = "", y = "Mean Deservingness Support", title = "")
ggsave("appendix_figure_a10_final.jpeg")








